  // Implements Firpo's estimator with fixed effects (not included in the pscore model)
  capture program drop qreg_col
  program define qreg_col
  {
   syntax varlist [if] , q(real) a(string)
   tokenize `varlist' 
   areg `1' `2' `3' `if', a(`a') 
   cap drop fe
   qui predict fe, d
   cap drop Ytem
   qui gen Ytemp=`1'-fe  
   ivqte Ytemp (`2') `if', continuous(`3') variance q(`q')
   drop fe Ytemp
  }
  end
  // Plot the results
  capture program drop quantplot
  program define quantplot
  {
  syntax, est(string) name(string)
  if "`name'"=="sound"    local titlename = "Knowledge of letter sounds"
  if "`name'"=="reading"  local titlename = "Fluency of oral reading"
  if "`name'"=="words"    local titlename = "Reading of non-words"
  if "`name'"=="readcomp" local titlename = "Reading comprehension"
  if "`name'"=="allliter" local titlename = "Literacy score"
  preserve
  keep if var=="`name'"
  #delimit ;
  twoway (rarea low up q, fcolor(gray%30) lcolor(gray%50) lpattern(dot)) 
       (line b q if est=="`est'" , lcolor(black) lwidth(medthick)) 
	   if est=="`est'", title(`titlename', color(black)) ytitle("") ylabel(0(0.2)1, 
	   noticks nogrid) xtitle("Quantile") xlabel(1(1)9) 
	   legend(off) graphregion(fcolor(white) lcolor(white) ifcolor(white) 
	   ilcolor(white)) plotregion(fcolor(white) lcolor(white) ifcolor(white) 
	   ilcolor(white));
  graph save Graph "$dir/results/plot_`est'_`name'.gph", replace;
  #delimit cr
  restore
  }
  end

  capture program drop quantileTE
  program define quantileTE, rclass
  {
  syntax varlist(max=1) [if], fixedeffects(string) treatment(string) controls(string) strata(string) [save(string)]
  // quantile treatment effects on sounds
   matrix qreg = J(9,3,0) 
   matrix firpo= J(9,3,0)
   local i = 1
   forvalues q = 0.1(0.1)0.9{
    //qte conditional
    qui qreg `varlist' `treatment'  `controls' `fixedeffects' `if', q(`q')
    matrix qreg[`i',1] = _b[`treatment']
    matrix qreg[`i',2] = _se[`treatment']
    matrix qreg[`i',3] = 2*(1-normal(abs(_b[`treatment']/_se[`treatment']))) 
    //firpo qte unconditional
    qui qreg_col `varlist' `treatment' `controls', a(`strata') q(`q')
     matrix firpo[`i',1] = e(b)
     matrix firpo[`i',2] = e(V)
     matrix firpo[`i',2] = sqrt(firpo[`i',2])
     matrix firpo[`i',3] = 2*(1-normal(abs(firpo[`i',1]/firpo[`i',2])))
	 local i = `i'+1
    }
	matrix Q = qreg \ J(1,3,.) \ firpo

	if "`save'"!=""{
	 preserve
     drop _all
     svmat Q
     drop Q3
     ren Q1 b
     ren Q2 s
     gen est = "qreg" if _n<10
     replace est = "firpo" if est==""
     drop if b==.
     gen q = _n
     replace q = q-9 if _n>=10
     gen low = b-(1.96*s)
     gen up  = b+(1.96*s)
     gen var = "`varlist'"
	 replace var = regexr(var, "_std_","")
     save "`save'", replace
	 restore
	}
	return matrix out = Q

}
  end
